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Abstract 

We use one-dimensional high-order central shock capturing numerical methods to study 
the response of various model solar atmospheres to forcing at the solar surface. The dy- 
namics of the atmosphere is modeled with the Euler equations in a variable-sized flux tube 
in the presence of gravity. We study dynamics of the atmosphere suggestive of spicule 
formation and coronal oscillations. These studies axe performed on observationally-derived 
model atmospheres above the quiet sun and above sunspots. 

To perform these simulations, we provide a new extension of existing second- and third- 
order shock-capturing methods to irregular grids. We also solve the problem of numerically 
maintaining initial hydrostatic balance via the introduction of new variables in the model 
equations and a careful initialization mechanism. 

We find several striking results: all model atmospheres respond to a single impulsive 
perturbation with several strong shock waves consistent with the rebound-shock model. 

These shock waves lift material and the transition region well into the initial corona, and 
the sensitivity of this lift to the initial impulse depends non-linearly on the details of the 
atmosphere model. We also reproduce an observed 3-minute coronal oscillation above 
sunspots compared to 5-minute oscillations above the quiet sun. 
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1 Introduction 

This paper investigates the response of the solar atmosphere to impulsive forcing at the so- 
lar surface using initial conditions derived from observationally-based atmosphere models. We 
simulate the response of model atmospheres above the normal solar surface and above sunspots. 
Various aspects of these simulations correspond to observed phenomena on the sun. Specifically, 
our simulations show behavior corresponding to spicules , linear features observed in the solar 
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atmosphere, and oscillations in the upper atmosphere. We find that the details of these phenom- 
ena depend on the atmosphere model in ways that correspond to the differences observed in the 


The dominant feature in our simulations is strong shock waves that propagate upward, lifting 
material from the lower atmosphere and causing oscillations of atmospheric materials. Because 
of these strong shocks, we perform our simulations using high-resolution central shock-capturing 
schemes. These schemes provide reliable approximations to solutions of the model equations in 
the presence of strong shocks while avoiding spurious numerical oscillations. We extend existing 
methods to computational meshes which have variable grid spacing. We also use new variables 
defined to facilitate numerical maintenance of hydrostatic equilibrium. 

One important result of our simulations is the observation that the height of the lifted material 
depends non-linearly on the details of the initial atmospheric model. A second result is that when 
we use the normal sun atmosphere model, the period of particle oscillations is in the five minute 
range while if we initialize our simulations using a sunspot model we see oscillations in the 
three-minute range. This strikingly reproduces observations of coronal oscillations in the solar 
atmosphere [2], 

Our simulations are based on the quasi-one-dimensional Euler equations applied to an initially 
hydrostatic atmosphere in a magnetic flux tube. We look at flux tubes that have constant area 
as well as a tube whose area increases with height. This model ignores heating, radiative energy 
loss and thermal diffusion, as well as magnetic fields. We find it remarkable that the naive model 
we use produces a quantitative match with observation. 

The structure of this paper is as follows: In section 1.1 we discuss basic properties of the 
solar atmosphere, describing its structure and summarizing the observed phenomena of interest 
to this paper. Section 2 presents the physical model that underlies our simulation. Section 
3 introduces our numerical method, including a discussion of initial and boundary conditions 
and the computational mesh. New second- and third-order reconstructions on irregular meshes 
are presented in section 3.2. In section 3.4 we present a new technique for maintaining initial 
hydrostatic equilibrium. Section 4 presents our results, focusing on the match between our 
simulations and observed properties of the solar atmosphere. 


1.1 The Solar atmosphere 

The solar atmosphere is a dynamic environment with high-energy phenomena occurring on many 
scales. At the base of the solar atmosphere is the photosphere , the surface of sun’s convective 
outer layer. The photosphere is roughly divided into two types of regions, quiet and active. The 
dominant features of the quiet photosphere are its temperature, about 5000 K, and granulations, 
which are currently understood as the surface of convective cells rising from below. The surface 
of the photosphere has been observed to rise and fall with velocities in the range of 1 km/sec, 
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Figure 1.1: Spicules near the solar limb, seen as darker linear features against the solar surface. The 
spicules appear darker in this narrow-band hydrogen-alpha image because their light is blue shifted due 
to their motion towards the observer. Photo credit: National Solar Observatory/Sacramento Peak 


and with a period of about 5 minutes. Active photosphere regions are associated with sunspots, 
which are regions of lower temperature, and very strong magnetic fields. Many complex and 
dynamic atmospheric features, such as prominences, magnetic loops and flares are associated 
with active regions. 

Above active regions magnetic flux tubes form, where a fixed number of magnetic field lines 
are confined by hydrodynamic pressure into a tube. Since hydrodynamic pressure drops with 
increasing height, the flux tube will in general expand with altitude, requiring a quasi- ID simula- 
tion. Though many flux tubes are curved and closed, following magnetic loops, we will consider 
the simpler case of open magnetic flux tubes pointing straight up. 

Above the photosphere, the solar atmosphere has a strongly stratified structure. The most 
striking feature of this stratification is the variation of temperature with height. For example 
above quiet regions, after slowly dropping in temperature to about 4700 K at about 3 5 meters, 
at 2 x 10 6 meters there is a strong rise to about 10 6 K within a few kilometers. This dramatic 
rise in temperature is called the transition region. The physical cause of the transition region 
is unknown. The solar atmosphere below the transition region is called the chromosphere, and 
above is the corona. 

It is well known that the vertically stratified structure of the solar atmosphere implies an 
acoustic cutoff-frequency that corresponds to a 5-minute period. Waves with periods shorter 
than 5 minutes cannot effectively propagate upward in the solar atmosphere. 
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Mdel Atmospheres 



Figure 1.2: Temperature profiles in the model atmospheres we use in our simulations 


One-dimensional models of the solar atmosphere have been developed, based on observed 
temperatures, densities and the requirement of hydrostatic equilibrium. The assumption of 
hydrostatic equilibrium is only approximately justified, since significant down- and up-flows have 
been observed. Models have been developed for the atmosphere above quiet regions [18] and 
sunspots [4, 11]. Fig. 1.2 shows the temperature profiles for the models that we use in our 
simulations. 

In the simulations described in this paper we begin with an atmosphere in approximate 
hydrostatic equilibrium. This initial atmosphere is based on one of the model temperature 
profiles described above. This atmosphere is then forced by an impulse perturbation at the base 
of the atmosphere. There are two observed phenomena in the solar atmosphere that correspond 
to features that appear in our simulations. 

• Spicules. Spicules [3, 14] are narrow near-linear features that are observed on the solar 
limb (see Fig. 1.1). They seem to be cooler photospheric material shooting up into the 
corona, and so can be modeled as an uplifting of the transition region in a magnetic flux 
tube. We will see such uplifting in our simulations. 

• Coronal Oscillations. Particle oscillations in the corona have been observed to have 
periods of about 5 minutes above quiet regions and 3 minutes above sunspots [2]. We 
observe a similar difference in oscillation period of coronal particles when we use the quiet 
sun as opposed to the sunspot model atmospheres. 
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We describe these observations in more detail in section 4 where we compare them directly 
with our simulations. 

1.2 Previous simulations 

There is a sizable literature of simulations similar to those presented in this paper, most of which 
attempt to model spicuie formation. Early attempts to model spicules as a ballistic uprise of 
the chromosphere into the corona due to a very large (60 km/sec) initial impulse are found in 
[13, 16]. These simulations used vertical straight flux tubes. Smaller impulses leading to a rise 
via the rebound shock phenomenon are considered in [1, 6, 7], some of which include expanding 
flux tubes. These simulations are extended in [15] by adding radiative energy loss, heating and 
heat conduction terms. 


2 The Model 

The quasi-one-dimensional hydrodynamic equations for material in a flux tube of area A (x, t ) in 
gravity are written in terms of density p(x,t ), velocity u(x,t ) and kinetic plus thermal energy 
E (x, t ) as 

(A P \ ( Apu \ / 0 \ 

Apu + Apu 2 = -Ap x - g (x) Ap + F (x, t) Ap , (2.1) 

\ AE J t \ A ( E +P) U J X \ ~g(x)Apu J 

where g (x) is the gravitational acceleration and the pressure p (x, t ) is given by the equation of 
state p = (7 — 1) (i? — \pu 2 ) with 7 = 5/3. The temperature T is given by p = R' pT where R! 
is the universal gas constant. F (x, t ) is a velocity forcing term which perturbs the base of the 
solar atmosphere in our numerical experiments. Various forcing terms will be considered. 

In order to use numerical methods that were developed for conservation laws, we rewrite (2.1) 
as a balance law (conservation law with source) 

( P \ ( pu \ ( -puA~ 1 A x - pA _1 A t \ 

( pu I + j pu 2 +p I = j -pu 2 A~' l A x — puA^At - g (x) p + F (x, t) p J. (2.2) 
\ E J t \C E + p)uJ x \ -(E + p)uA- l A x -EA~ l A t -g{x)pu j 

This form makes explicit that the scale of the flux tube area does not effect the dynamics: only 
the change in area enters the equation. Following previous authors [6, 15] we consider static flux 
tubes and have omitted the the A t terms from (2.2). 
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For our flux tube geometry we adapt the area function of [15]. We use the functional approx- 
imation 

A ( x ) = Aq [1 + 6 (a ( x ) — Oo)] 2 , (2-3) 

where b = 3.2, a(x ) = tanh (s (x — xq)) with s = 10 -6 and xq = 1.5 x 10 6 , and ao = a(0). 


Aq — A (0) is the area of the flux tube at the base of the atmosphere, but does not appear in the 


simulation. Therefore 

A' (s) _ Zbs (1 - a (x) 2 ) 
A(x) 1 + b (a ( x ) — ao) ' 


3 The Numerical Method 

In this paper we use high-order, non-oscillatory central schemes designed to solve systems of 
conservation laws with source terms of the. form 

Qt + f(q) x = S(q,x) . (3.1) 

Here q G R p is a p-dimensional solution vector, / is the flux function, and S' <G K p is the p- 
dimensional vector of source terms. The solution of 3.1 may become singular in finite time, 
which in turn requires careful study when dealing with numerical approximations. 

Our approach for solving (3.1) is based on extending the semi-discrete central-upwind scheme 
of Kuragnov et. al. (KNP) [9] to irregular grids. The KNP method is a simple efficient scheme 
with less dissipation than fully central schemes, and has desirable stability properties. Any such 
scheme is composed of a numerical flux, a piecewise-polynomial reconstruction and an ODE 
solver. The order of accuracy of the reconstruction and the ODE solver determines the order of 
accuracy of the method. Below we address all three ingredients, with a particular focus on our 
new high-order reconstructions on irregular grids. 

Throughout this section we assume a one-dimensional grid with irregularly spaced nodes 
Xj. The distance between consecutive grids points is denoted as A Xj := Xj + 1 — Xj. We define 

x j+ i := (xj + Xj + i) /2 and the cell Ij = j . For any function / (x) we use the notation 

fj := f (xj). The cell average of q in the cell Ij is given by 
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We assume that the cell-averages are known at time t n . The first step in the derivation 
of the approximate solution is to generate a piecewise-polynomial reconstruction from these cell- 
averages. Such a global reconstruction is defined as 

H x ) = J^q j (x)xi j (x), (3.2) 

j 

where xi } ( x ) is the characteristic function of Ij, and <jj (x) are polynomials of a suitable degree. 
The construction of admissible qj will be discussed in the section 3 . 2 . For the time being, it is 
sufficient to assume that a reconstruction of the form (3.2) is known. 

We denote the point- values of q at the interfaces of the cell Ij by 



3.1 The KNP central-upwind method 

It is straightforward to extend the semi-discrete central scheme of [9] to irregular grids. When 
written in a conservative form this scheme takes the form 


dqj H :+\ ~ H ]-\ = 

dt x j+ i — Xj_i j ’ 

J-t- 2 J 2 


(3.3) 


where Sj is a discretization of the averaged value of the source term in Ij and the numerical flux 
# 7 + 1/2 is given by 


H j+i 


<\ f ( q 7+h) + a 7n f ( q U) 


+ Q ^ 

a iH +a J+l 


a + , , a . , 

3+2 3+ 


a + ! 

3+1 


3+7 


The local speeds of propagation of information from the discontinuities on the interfaces of the 
cells, a,j± i, are given by 


a + , i — max 

3+2 






min ( A ‘(f)’ 0 ) ■ 


Here \k denote the eigenvalues of the Jacobian of / evaluated at Xj + 1 . In the case of 
problem ( 2 . 2 ), we have 


at ,i = max(u — c,u,u + c, 0 ) , a. . 


| min (u — c,u,u + c, 0 )| , 
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where c = \rypjp is the speed of sound and all fields are evaluated at ar J+ i . 

We approximate the cell- average of the source term Sj as S (pj,Uj, Ej,Xj). This amounts 
to a first-order quadrature approximation to the cell average. The O (A Xj) error from this 
approximation appears in the terms containing the flux tube area variation and the gravitational 
acceleration g (x). This error is acceptable since the flux tube model (2.3) is itself only qualitative. 
While we compute the gravitational acceleration g (x) for each cell, it is slowly varying and so is 
approximately constant in each ceil. 

3.2 Reconstructions and order of accuracy 

In this section we introduce new second- and third-order reconstructions for irregular grids that 
are used in the scheme (3.3). In spite of the uncertainties that exist in the model atmospheres 
we use for our initial conditions as well as the overall simplified nature of our model, high-order 
methods are valuable since they enable us to conduct long-time simulations with well resolved 
shock structures. 

Our second-order reconstruction, when used in (3.3), is highly efficient and leads to a method 
that is total variation diminishing (TVD). The TVD property assures that spurious numerical 
oscillations will not increase from timestep to timestep. The third-order method, when used in 
(3.3) provides higher resolution and leads to a ’’number of extrema diminishing” (NED) method 
on regular grids. While the NED property does not guarantee that existing oscillations will not 
grow, it does assure us that new numerical oscillations will not be introduced. Though we have 
not demonstrated that the NED property holds for irregular grids, the computational mesh we 
use in our simulations is regular in the region of greatest interest. 

The reconstruction qj(x) in the cell Ij should satisfy the following requirements: 

• Conservation. 



While any linear reconstruction will have this property on a regular grid, this will not be 
automatically true when the grid has variable spacing. 

• Accuracy, qj is an rth-order accurate reconstruction, i.e. if we scale A Xj — ► SAxj then 
the error scales as S r . On a regular grid with grid spacing Ax this amounts to 

Qj ( x ) = q (x) + O (Ax r ) , x e Ij. 

• Oscillation minimizing. Spurious numerical oscillations should be avoided. 
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3.2.1 A second-order reconstruction 


Second-order reconstructions are provided by piecewise-linear interpolants, with a different in- 
terpolant for each cell Ij. A linear interpolant that is conservative in Ij for any Dj is given 
by 

a. (x\ = <7„- -I- D.- ( , x — ( X; -4- — fA.r.- — . (3.4) 

v ' > ■> y y ■> 4 ' J J *' j j 

The discrete derivative, Dj, is taken as 

Dj = MinMod ( 4- , 4- ) . (3.5) 

\ Axj + 2A Xj—i + Axj — 2 Axj+i + 2Axj + Axj_i J 

where the MinMod limiter is defined as 

{ min (?i,..., 5 „), $ > 0,Vz, 

max (?!,..., q n ), q { < 0, Vi, 

0, otherwise. 

See [17] for a discussion of nonlinear limiters on uniform grids. The expression (3.4) now provides 
the required values at the cell interfaces: 


4,+ ! = <h ) =9j + (A Xj + Axj-J , 

Qj + 1 = %+i (^+i) = ty+i - (&cj + Ai j+ i) . 

Lemma 3.1 The reconstruction (3.4) is total variation preserving. 

Proof: Our reconstruction is R(v) = V. qj(x)xj(x) : where qj is given by (3.4). Then 


TV(R(v)) = V ( Ax ' + Ali - 1 


+ <7j - 9j+i + ^(A^j+i + Axj) H — —^(A Xj + Axj+i) 

It is easy to verify that Dj defined in (3.5) satisfies 

1. sgn (Dj) = sgn(D J+1 ) = sgn($ +1 - %) 

2. | KAxj^+A Xj)Dj + (A xj + Ax j+1 )D j+1 \ < \q j+1 - qj\. 
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Then (assuming periodic boundary conditions) 

, ( A Xj + A 

- L'^'I V 2 

3 L x 


TV(R(v)) = 


+ I Qj - qj+il - ^(Aajj+ 1 + Axj) - ^^(Axj + Axj+i) 


— y~! I Qj Qj+ 1 1 T l-Dj 


z\x,_i + zix 




= Tl/(u). 


3.2.2 A third-order reconstruction 

The third-order method from [8], based on [12], uses a linear combination of a linear interpolant 
Lj (x) and a quadratic interpolant Qj (x) that is third-order accurate in smooth regions. When- 
ever a. discontinuity is present, the order of the reconstruction is reduced by choosing one linear 
interpolant from the appropriate side of the discontinuity. This latter choice minimizes spurious 
numerical oscillations. The reconstruction of [8] is given in terms of the convex combination 

Qj (x) ■■= (1 - Oj) Lj (x) + 6jQj (x) , (3.6) 


where 


Oj ~ < 


M 


mm 


mm 






j 


M 


l i, 


-r L K x i-h) ? 

fj-Lj (V,_^ j mj-L] (xj + ^ 


m i+r L A^±i 


Qj - 1 ^ Qj < Qj+ 1) 


, lj , 5j_i > qj > q j+ 1, 
otherwise, 


and 


Mj = max {Qj (x J+ i) , Qj (xj_Aj ] , m ] = min { Qj (x j+ i) , Qj (xj_ i) } , 
M i±i = m ax ( Lj (x J± i) + L j±1 (^±|)) , Qj ± i (x,±i) | , 


m j± i_ = min ( Lj (x j± Aj + L j±1 (x J± i)) ,Q j± i (^±i)} • 


12 


S. Bryson, A. Kosovichev and D. Levy 


as 

by 


We now extend the method (3.6) to non-uniform grids. The linear function, Lj (x), is taken 
the conservative reconstruction (3.4) with (3.5). The quadratic function, Qj(x), is replaced 


Qj (x) = Aj + Bj (x - (xj -I- x 0 )) + Cj (x - (xj + x 0 )) 2 , 


(3.7) 


where 

x 0 = ^ (A Xj - Axj-i ) , 

Aj = Qj ~ ^g7 jCjt 

Bj — 2 otj [ f3j (Axj_i T 3 A Xj + 2Axj+i) Qj— i T PjPj+iftjQj 
+ Pj+i (2Axj-2 + 3Axj_i + A Xj) Qj+i ] , 

Cj — 12a j [(3jQj-i + Pjpj+irjjQj + Z^+i^+i] • 

Here jj = A Xj + A x 3 _x, 

aj (Axj-2 + 2Axj-i + 2Axj + Axj + i) ’ Axj-2 + 2Axj_i + Axj ’ 

Tjj — — (Axj_2 + 3Axj_i + 3A Xj + Axj^-i) , 


y.j = (2Axj_ 2 -I- 5Axj_i -)- 5A Xj + 2A j+ i) (—Axj-i — Ax^-i + Axj + Ax J+ i) . 

It is easy to verity that the interpolant (3.7) is conservative and third-order accurate on I r 
Using (3.7) we get that the point values at the cell interfaces are 

Qj (x j+ 1 ) = -ajpjjjij+iqj-! + ^1 + latjpjPj-ijj (TjW + Mj)') Q + a j^p~^+ 


Qj (^-i) = + { 1 + \ a jP}0i-rt3 (Wi ~ &)) q j ~ +i7i7j- 1 9i+i • 


and 
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3.3 Time integration 

The time stepping is performed by either a second- or third-order TVD Runge-Kutta method 
from [5]. Given data q T - at time t n , we advance to the next time step by solving an ODE of the 
form 

-% -Gf,). (3-8) 

A second-order TVD-RK method for (3.8) is 
q (1) = q n + AtF(q n ), 

?” +1 = ilV + ^’ + AtF ( g W)). 

A third-order TVD-RK method for (3.8) is 
<7 (1) = q n + AtF(q n ), 

9® = \< r + 1 («<■>+ (,(«)) , 

9" +1 = i«” + ! (9® +AiF (,<»)). 

Each time step is set to 

At = /3max 

j a (Jj) 

where a (Jj) is the spectral radius of the Jacobian J = df/dq evaluated at Xj. The parameter f3 
is taken in our simulations as P = 0.9. 


3.4 Hydrostatic Equilibrium 


We assume hydrostatic equilibrium at initialization and at the computational domain boundaries. 
This implies that u = 0 and there is no time evolution so long as hydrostatic equilibrium is 
maintained. Therefore hydrostatic equilibrium relates pressure and density via (2.2) as 


dp ( x ) 

dx 


= ~9 (x)p(x) = - 


ff (z)p(t) 
R'T(x) ‘ 


(3.9) 


We use the temperature profile T ( x ) and pressure at the bottom of the atmosphere p (0) = 
R'p (0) T (0) that are given by an atmosphere model to solve (3.9) for the hydrostatic pressure 
profile p (x). 
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3.4.1 Hydrostatic Balance 

In order to obtain numerical hydrostatic balance, we modify the model (2.2) through a change of 
variables. We assume that we are given an initial density po i x ) and energy Eq (x) in hydrostatic 
equilibrium as described in section 3.4.2. We define the variables p := p — po and E := E — Eo 
and write the system (2.2) in the mathematically equivalent form 


(p + Po) u 
(p + po) u 2 + p 


E + Eq 


( i „ \ 4-1 A 

~r HO) usx jx x 


— {p + Po) (u 2 A l A x + g ( x ) — F (x, t )) 

- (e + Eq + pj uA~ l A x - g (x) (p + p 0 ) u 


, (3.10) 


with p (7 — 1) ^ E + E 0 — | (p + p 0 ) u 2 ^ . Here we define u := pu/ (p + po) and use the fact 
that hydrostatic equilibrium implies that the flux tube is static so A t = 0. With these new 
variables, p — E = 0 as long as the system remains in hydrostatic equilibrium, i.e. satisfies (3.9) 

with u = 0. We estimate po and Eq at the cell interfaces using the interpolants 

described in section 3.2. As described in section 3.4.2, if po (x) is defined by (3.12), then our 
numerical method (2.2) will remain in hydrostatic balance as long as u = 0. 


3.4.2 Hydrostatic Initialization 

We use (3.9) to initialize our simulations in hydrostatic equilibrium given a temperature profile 
specified by a Solar model. Equation (3.9) is integrated on a high-resolution grid using the stan- 
dard fourth-order Runge-Kutta method, then linearly interpolated onto the actual computational 
grid of the simulation. This process gives an initial pressure profile po (x) which is transformed 
into an energy profile via E 0 (x) = po(x) / (7 — 1). 

The initial density profile po (x) is determined by the requirement of numerical hydrostatic 
equilibrium u = 0 in our method (3.3). In this case we have of ' 1 = a 7 , 1 = c.+i. Inserting these 

3~r 2 J~T~ 2 J 2 

conditions into (3.3) applied to our model (2.2) (supressing the subscript 0 on p 0 for clarity), we 
obtain 


2 x i- 


P ~P 

p~ +p + 

E~ — E + 


( P ~P 

p~+p+ 

\ E~-E+ 


(3.11) 


We therefore have hydrostatic balance in the second component if we initialize the density term 
in the source as 


(Po)j = 


2 9i K+i - x j- 


EH +p ti-( p J-i +p U). 


(3.12) 
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The first and third components of (3.11) will vanish if we use the variables defined in section 
3.4.1 so long as hydrostatic equilibrium is satisfied at Xj. Therefore if we initialize p according 
to (3.12), p = pu = E = 0 and Eq — Po/ (j — 1), then the method (3.3) applied to (3.10) will 
maintain hydrostatic equilibrium. 


Remark. Eq. 3.12 is a natural discretization of the cell average using the cell interface values 

Pj±V 



which follows from the condition for hydrostatic equilibrium. Eq. 3.12 cannot be used to ap- 
proximate the cell average of the source term after initialization because the system does not 
remain in hydrostatic equilibrium after we apply our forcing term. 


3.5 Boundary conditions 

Even though we use stretched grids, placing our boundaries at great distances as described in 
section 3.6, we want to limit that distance so we can use as few grid nodes as possible in the 
computational domain. We therefore have two issues to consider: the failure of initial hydrostatic 
equilibrium at the boundaries, and the problem of reflections, particularly of strong shocks. 
Failure of initial hydrostatic equilibrium at the boundaries will cause a wave to propagate into 
the domain. 

Zeroth-order boundary conditions, where qj = q\ for j < 1 and q 3 = <jv for j > N where N is 
the number of grid nodes, prevent reflections off the boundary but violate hydrostatic equilibrium. 
We therefore adapt a dynamic boundary condition approach, where the boundary conditions are 
hydrostatic until a high velocity value propagates from the interior to that boundary: if u has not 
risen above a threshold U since the start of the simulation then hydrostatic boundary conditions 
are applied. Otherwise zeroth-order boundary conditions are applied. This approach allows small 
reflections, which propagate more slowly than strong shocks thereby maximizing the time until 
the reflections enter the lower atmosphere region. For a typical simulation U is taken to be 1 
km/sec. 

Hydrostatic boundary conditions are implemented by defining the field values in the ghost cells 
via linear interpolation of the high-resolution initialization grid described in section 3.4.2. This 
requires that the high-resolution initialization grid’s domain be larger than the computational 
domain plus any ghost cells. 
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3.6 The computational domain 

The computational domain is divided into three regions: the lower atmosphere, which is our 
domain of primary interest and represents the chromosphere and lower corona, the upper corona, 
and a non-physical region below the solar surface. The lower atmosphere region has constant 
high-density grid spacing while the upper atmosphere and sub-surface regions have exponentially 
stretched meshes. The stretching of each stretched region is defined so that grid spacing matches 
that of the (non-stretched) lower atmosphere at the shared boundary, and extends to a specified 
distance with a specified number of nodes. In this way we can have high resolution in the region 
of interest and boundaries sufficiently far away such that reflected waves do not enter the lower 
atmosphere region during the course of the simulation. Typical simulation values place the upper 
boundary of the lower atmosphere at 10 7 meters and the upper corona region at 5 x 10 8 meters. 

The lower atmosphere domain is defined so that the solar surface is located at a coordinate 
value of 7.5 x 10 6 meters, which corresponds to a solar radius of r© = 6.9599 x 10 8 meters. Thus 
x = 0 corresponds to r 0 := 6. 8849 x 10 8 meters. 


The non-physical interior region. In the non-physical interior region the physics is adjusted 
to delay reflections from the boundary so they do not enter the physical region during the course of 
the simulation. We define the gravitational acceleration in the interior region to be exponentially 
damped with increasing depth: 

1, x > 0, 

ex p{-(a^) 2 }> * <0 > 


9 O) = 


GMq4> (x) 


<P (*) = s 


where M© is the solar mass, r = r Q + x and xj is a scale defined so that g (x) ~ 0 in most of the 
interior region. In a typical simulation Xj is typically taken as — 1 x 10 7 meters. 

In hydrostatic equilibrium we have dp/dx = — gp = —gpj ( R'T ) so a smaller g will lead to a 
smaller growth in p. We extend the temperature T into the interior as a constant: T (x) — T (0) 
for x < 0. The result is a slower speed of sound so the bottom of the interior region can be closer, 
requiring fewer grid nodes, while still avoiding reflections off the lower boundary from entering 
the lower atmosphere region. 


4 Results 

In this section we present various results using the model and methods in this paper. We are 
particularly interested in the response of the solar atmosphere to a single impulse at the solar 
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surface. Our focus is on simulations that bring out phenomena that resemble those observed in 
the solar atmosphere. 

We consider the forcing term F (x,t) in (2.2), defined as 
F(x,t) = ^^X(x)T(t), 

where A is a constant velocity amplitude in m/sec. X and T are square functions with respective 
widths w x and Wt- for |t — to\ < w t , T (t) := 1, otherwise 0, and for \x — Xo\ < w x , X (x) 1, 

otherwise 0. 

In the above, A,w t and w x are tunable parameters that together determine the total energy, 
deposited into the atmosphere by the forcing function. We will see below that simulation results 
can depend strongly on these parameters, though we find phenomena matching those observed 
for physically reasonable values of these parameters. Our general approach will be to tune the 
parameters so that some aspect of the simulated phenomena matches observations and compare 
other aspects of the simulation. For example, we may choose parameters that cause the transition 
region to rise to around 5 x 10 6 meters and study the resulting particle oscillations, atmospheric 
profiles and shock wave development. For simplicity we assume that the solar atmosphere is 
entirely composed of hydrogen. 

We will see below that the specific shape of the initial atmosphere may have a strong impact 
on the resulting phenomena. 

4.1 Accuracy and comparison of the reconstructions 

We first check the accuracy of the second- and third- order interpolants described in section 

3.2 on an irregular grid. Table 4.1 shows the relative L 1 and ''relative L°° errors for the linear 
advection problem 

&+9rr = 0 (4.1) 

with periodic boundary conditions. The errors are shown at T = 1 for the periodic initial data 
q (x, t = 0) = sin 4 (wx) on the domain [—1, 1] with nodes at Xj = Xj- 1 + A (1 + sin (2wj/N) /l.l), 
where N is the number of grid nodes. Here A = 1/(7V — 1). 

We now compare the quality of the second- and third-order reconstructions for our simulations 
of the solar atmosphere. Fig. 4.1 shows the result at T = 570 seconds, initialized using the quiet 
sun model and a single-timestep square impulse at T = 0. Two simulations, with the second- 
and the third-order methods, were performed with a low resolution of 200 grid nodes in the 
lower atmosphere region in order to highlight the difference between the second- and third- 
order methods. The other two simulations were performed with 1000 grid nodes in the lower 
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Second- order Reconstruction 

Third- order Reconstruction 

N 

L L Error 

Rate 

L°° Error 

Rate 

L 1 Error 

Rate 

L°° Error 

Rate 

50 

0.332 

- 

1.97 x 10“* 

- 

8.69 x 10"-* 

- 

5.33 x 10" 3 

- 

100 

0.100 

1.73 

4.72 x 10~ 3 

2.06 

1.37 x 10“* 

2.66 

4.17 x 10~ 4 

3.68 

200 

2.95 x 10"* 

1.76 

1.02 x 10“ 3 

2.20 

2.20 x 10” 3 

2.64 

3.26 x 10" 5 

3.68 

400 

8.93 x 10- 3 

1.72 

2.18 x 10~ 4 

2.23 

3.79 x 10" 4 

2.54 

3.65 X 10~ 6 

3.16 


fable 4.1: Relative L l and errors and convergence rates for the reconstructions on irregular 
grids used in the scheme (3.3) applied to the advection equation (4.1). 




Figure 4.1: Comparison of second and third order methods. Left: Velocity. Right: Temperature. 

atmosphere region. We see that there is a significant difference between the two methods in the 
low-resolution case, with a much smaller difference in the high resolution case. Also notable in 
this figure is the strength of the shock in velocity value. It is these large velocity shocks that 
cause the uplift and particle oscillations. 

Based on these results we will perform the rest of the simulations in this section with 1000 
grid nodes in the lower atmosphere region using the second order method. 

4.2 Spicule formation 

Spicules are narrow near-linear features that are observed on the solar limb [3, 14] (see Fig. 1.1). 
They are observed to rise with upward velocities of about 25 km/sec, have temperatures in the 
range of 5000 to 10,000 K and densities of about 3 x 10 -13 g/cm 3 . Their maximum height 
is between 6.5 x 10 6 and 9.5 x 10 6 m. Because of these relatively low temperatures and high 
densities, spicules can be thought of as an intrusion of chromospheric material into the corona. 
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Figure 4.2: History of the temperature {left) and density (right) profile for the quiet sun model in 
response to a square, single-timestep 1000 m/sec impulse. 

Their linear shape can be explained as this material being channeled inside a magnetic flux tube. 

We are primarily interested in the behavior of the transition region in response to impulses 
at the base of the atmosphere, as well as the physical conditions below the transition region 
at its time of maximum height. We define the position of the transition region to be where 
the temperature rises through 200,000 K. We compare simulations using the three atmosphere 
models with both constant and expanding flux tubes. We do this despite the fact that spicules 
are not observed above sunspots in order to study the dependence of spicule formation on the 
atmosphere model. 

The values of w t and w x in the impulse function F (x, t ) (see the introduction to section 4) 
are chosen so that for A = 1000 m/sec the transition region is lifted to just above 5 x 10 6 meters 
in the quiet sun atmosphere. These values of w are different for the case of constant flux tube 
width vs. an expanding flux tube. For a constant flux tube we use w t = 0.065 seconds and 

= 1.4 x 10 5 meters. For expanding flux tubes significantly longer impulses were required to 
generate the required lift of the transition region: in this case w t = 0.325 seconds. 

The dominant feature of our simulations is strong initial and rebound shock waves that re- 
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peatedly appear for as long as 30 minutes after the single initial impulse. An example of the 
effect of these shocks on the quiet sun model is shown in Fig. 4.2. The shock trajectories for 
A = 1000 m/sec are shown for the various model conditions in Fig. 4.3. We see a character- 
istic increase in speed as the shock crosses the transition region from the chromosphere to the 
corona. As time passes, this increase in speed occurs at higher altitudes as the transition region 
rises. The expanding flux tube also shows short-lived downward propagating shocks below the 
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by the time they reach the transition region. This causes a dramatic uplift, i.e. lifting mate- 
rial and the transition region well into the initial corona. We track the motion of atmospheric 
material by seeding test particles at specific heights, then tracking their motions throughout the 
simulation. The position x (£) of these particles is updated at each time step by solving the ODE 
x f (£) — u (x, t ) via the standard fourth-order Runge-Kutta method, where u (rr, t) is the velocity 
component of the approximated solution. The resulting particle histories is shown in Fig. 4.4. 

We see immediately from Fig. 4.4 that the Lites-Skumanich sunspot model is much more 
responsive to the initial perturbation, with much higher uplift than the other two models. This 
is striking since, as seen in Fig. 1.2, the Lite-Skumanich model is in some sense intermediate 
between the other two models. This demonstrates a sensitive non-linear dependence of the ability 
to form spicules on the details of the atmosphere model. This behavior is also evident in Fig. 4.5, 
which show's the maximum height of the transition region for varying impulse amplitudes. The 
three models exhibit a near-linear dependence on impulse amplitude for the constant flux tube, 
while for the expanding flux tube the Lites-Skumanich sunspot model shows a greater sensitivity. 
The Ding-Fang sunspot model is least responsive, requiring stronger impulses to rise to desired 
heights,. We see that it is not difficult to raise the transition region to 5 x 10 6 meters, the height 
appropriate for spicules, in all cases, though we remind the reader that we tuned our parameters 
for such a rise in the quiet sun model. 

Fig. 4.6 shows estimates of the speed of the transition region in its initial rise for varying 
impulse amplitudes. This estimate is defined as the height of the first local maximum divided by 
the time taken to reach that height from first motion. In all cases the speed of the rise increases 
approximately linearly with amplitude. For the constant flux tube with impulse amplitudes of 
1000 m/sec we see rise speeds of about 22 km/sec for the quiet sun and Lites-Skumanich sunspot 
models, while the Ding-Fang sunspot model shows a slower rise. For an expanding flux tube we 
see somewhat slower rise speeds, about 18 km/sec for the quiet sun and Lites-Skumanich models 
and 1000 m/sec impulse. 

To examine the physical conditions below the transition region at time of maximum rise we 
consider impulse amplitudes that lift the transition region just above 5 x 10 6 meters, which is 
different for each atmosphere model. The temperature and density profiles for these cases are 
shown for the three atmosphere models and both flux tube types in Fig. ??. The properties 
observed from this figure are summarized in table 4.2, where we see that our simulations match 
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Figure 4.3:- The trajectories of the velocity shocks for an impulse amplitude of 1000 m/sec. Left: 
Constant flux tube. Right: Expanding flux tube. Top: Quiet sun model. Middle: Ding-Fang sunspot 
model. Bottom: Lites-Skumanich sunspot model. Note the downward propagating shocks in the 
expanding flux tube. 
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Figure 4.5: Maximum height of the transition region for various surface velocity perturbation ampli- 
tudes. Left: Constant flux tube. Right: Expanding flux tube. 



Figure 4.6: Average speed of rise of the transition region from initial position to first peak height vs. 
impulse amplitude. Left: Constant flux tube. Right: Expanding flux tube. 
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the observations remarkably well. 



Constant flux tube 

Expanding flux tube 


Obseved 

VAL 

DF 

LS 

VAL 

DF 

LS 

Temperature (K) 

0.5 x 10 4 - 10 4 

2 x 10 4 

2 x 10 4 

2 x 10 4 

10 4 

10 4 

3 x 10 4 

Density (g/cm 6 ) 

3 x 10~ 13 

2 x 10" 13 

io ~ 14 

10“ 14 

3 x 10~ 13 

2 x 10" 14 

10" 14 


Table 4.2: Temperature arid density below, the transition region 
on Fig. 4.7 compared with observed spicule properties. 
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4.3 Coronal oscillations 

Recent observations in [2] have detected a notable difference in the period of oscillation of coronal 
particles above sunspots compared to periods above the quiet sun. Above sunspots, the observed 
oscillations have periods clustered around 3 minutes, while above the quiet sun the periods are 
clustered near 5 minutes (see Fig. 4.8). 

In our simulations, we find a similar relationship between period of oscillation of particles 
in the corona (above the transition region). As seen in Fig. 4.4, coronal material exhibits an 
oscillatory behavior after being lifted by the shocks. Due to the multiple shocks this behavior is 
complex, but there are dominant modes of oscillation. 

For all model atmospheres, the period of oscillations of our coronal particles depends on the 
amplitude of the initial impulse, though the details of this dependence differs from model to 
model. We study of this dependence by fixing the impulse parameters w t and w x for each model 
and varying A. We use the same values for w t and w x as in our spicule study. 

The coronal oscillations found in our simulations are summarized in Fig. 4.9 and discussed in 
detail below. Example power spectra for A = 1000 m/sec are shown in Fig. 4.10. We find that 
there is a significant difference between the case of a constant width flux tube and an expanding 
flux tube. 

4.3.1 Constant flux tube 

In Fig. 4.9 we see that when we use the quiet sun model atmosphere as our initial condition, 
we find a monotonically increasing relationship between the initial impulse amplitude and the 
period of coronal oscillations. The period of peak oscillation in this case ranges from 4 minutes 
for a 200 m/sec impulse to 6.5 minutes for a 1200 m/sec impulse. 

When we use the Ding-Fang model atmosphere above sunspots, we find a more complex 
situation. For low-amplitude 200 to 400 m/sec impulses, we find coronal oscillations with periods 
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Figure 4.7: The temperature (solid line) and density (dashed line) profiles at the time of maximum 
transition region height in response to a square, single-timestep impulse. Left: Constant flux tube. 
Right: Expanding flux tube. Top: Quiet sun model. Middle: Ding-Fang sunspot model. Bottom: 
Lites-Skumanich sunspot model. 
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above sunspot 



oscillation period (minutes) 


Figure 4.8: Histograms of observed periods of coronal oscillations. Top: Above the quiet sun. 
Bottom: Above sunspots. 



Figure 4.9: The simulated periods of oscillations for the three atmosphere models and various surface 
velocity perturbation amplitudes. Left: Constant flux tube. Right: Expanding flux tube. 
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in the 3.5 to 5 minute range. For an amplitude of 700 m/sec, we find a bimodal regime with one 
mode of oscillation at about 3.2 minutes, and the other at about 5 minutes. For amplitudes from 
900 to 1000 m/ sec the longer-period mode is significantly damped, and the dominant mode has 


periods around 3.5 minutes. 

The Lites-Skumanich sunspot model shows a monotonic dependence on impulse amplitude, 
but the dominant mode of oscillation jumps from between 3 and 3.5 minutes, for amplitudes 
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oscillation increases almost linearly for impulse amplitudes from 500 to 1200 m/sec. 

We see that for a constant width flux tube the Ding-Fang sunspot model and the quiet sun 
model shows coronal oscillations that are consistent with observed coronal oscillations above 
sunspots and above the quiet sun. 


4.3.2 Expanding flux tube 

For expanding flux tubes we see in Fig. 4.9 a bimodal behavior for all the atmosphere models. 
In the quiet sun model we find the dominant mode of oscillation increasing in the four minute 
range for impulse amplitudes ranging from 300 to 600 m/sec. For impulses between 700 and 1000 
m/sec we see bimodal behavior, with the longer mode’s period rising to over 5 minutes, and a 
shorter mode in with a near-constant period of 3.5 minutes. For higher amplitudes the slower 
mode is damped and the dominant mode is just over 3.5 minutes. 

In the Ding-Fang sunspot model we have a single dominant mode for impulse amplitudes 
between 300 and 1000 m/sec, with the period rising monotonically from 3.5 to just under 5 
minutes. At an amplitude of 1100 m/sec we see bimodal behavior, with a slow mode of 5 
minutes and a fast mode of 3.5 minutes. At 1200 m/sec the slow mode is damped and the 
dominant mode is about 3.5 minutes. 

For the Lites-Skumanich sunspot model, we see a bimodal behavior at the lower amplitude of 
300 m/sec, with 3 and 3.8 minute modes, and amplitudes of 1000 m/sec and above. In between 
we have a dominant mode which increases from about 3.3 minutes to 4 minutes. For amplitudes 
of 1000 m/sec and above we have a very slow mode of over 8 minutes and a fast mode of just 
over 4 minutes. 

4.4 Effect of Moving the Transition Region 

In an attempt to understand the role the location of the transition region plays in the above 
results we present simulations in which the quiet sun model is scaled so that the transition region 
has various heights. We use the standard square 1000 m/sec impulse for this study. 

In Fig. 4.11 we see show the dependence of the maximum transition region height on the 
initial transition region height. We see that there is a non-linear dependence, where initially 
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Figure 4.10: Power spectra of particle oscillations for a square, single-timestep impulse. The spectra 
are stacked in order of the height of the particles. Left: Constant flux tube. Right: Expanding flux 
tube. Top: Quiet sun model. Middle: Ding-Fang sunspot model. Bottom: Lites-Skumanich sunspot 
model. 
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transition region maximum height, constant flux tube 



Figure 4.11: Maximum height of the transition region in response to a square impulse of 1000 m/sec 
as a function of the initial transition region height in a scaling of the quiet sun model. 

higher transition regions get lifted much higher. One possible explanation is the fact that the 
rising shocks increase in strength as the density decreases in the hydrostatic atmosphere, which 
results in higher velocities by the time the shock reaches the transition region. 


5 Conclusions 

In this paper we used a high-order shock-capturing scheme to study the response of various 
one-dimensional model solar atmospheres to a single impulse at the solar surface. Our treat- 
ment includes an extension of existing methods to irregular grids and a new technique which 
numerically maintains initial hydrostatic balance. In spite of the simplicity of our model we 
find a remarkable correspondence between our results and observation: 3-minute coronal oscil- 
lations above sunspots compared with 5-minute periods above quiet regions observed in [2] and 
appropriate spicule properties [3, 14]. These results suggest that a viable mechanism of spicule 
formation is large velocity shocks in response to surface impulses. We further find that the lift- 
ing of material in response to an initial impulse depends non-linearly on the details of the initial 
atmosphere model. 

The next steps in our study include adding energy loss and heating to our model follow- 
ing [15]. This will allow us to study the atmospheric response to continuous surface forcing. 
Extensions to two dimensional studies will follow, where the model will be extended to the 
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magnetohydrodynamic equations. 

Acknowledgment: The work of D. Levy was supported in part by the National Science Foun- 
dation under Career Grant No. DMS-0133511. 


R ofidron fnc 

V' ' JL. > ' JU 

[1] Andreev, A., S., Kosovichev, A. G., On a Mechanism of Spicule Formation by Shock Waves 
in Magnetic Tubes , Astronomy Letters, 403, (1994), pp.323-326. 

[2] De Moortel, I., Ireland, J., Hood, A. W., Walsh, R. W., The Detection of 3 & 5 Minute 
Period Oscillations in Coronal Loops, Astronomy and Astrophysics, 387 , (2002), pp.L13- 
L16. 

[3] de Pontieu, B., PhD Thesis, (1996), University of Ghent, Belgium. 

{4] Ding, M. D., Fang, C., A Semi-Empencal Model of Sunspot Penumbra , Astronomy and 
Astrophysics, 235 , (1989), pp.204-212. 

[5] Gottlieb S., Shu C.-W., Tadmor E., Strong stability-preserving high order time discretization 
methods , SIAM Review, 43 (2001), pp.89-112. 

[6] Hollweg, J. V., On the Origin of Solar Spicules, ApJ, 25 7 , (1982), pp. 345-353. 

[7] Kosovichev, A. G., Popov, Yu. P., The Computation of One- Dimensional Non- Stationary 
Problems of Gravitational Gas Dynamics, USSR Comput. Maths. Math. Phys., 19, (1998), 
pp. 168-175. 

[8] Kurganov, K., Petrova, G., A Third-Order Semi-Discrete Genuinely Multidimensional Cen- 
tral Scheme for Hyperbolic Conservation Laws and Related Problems , Numerische Mathe- 
matik, 88, (2001), pp. 683-729. 

[9] Kurganov, K., Noelle, S., Petrova, G., Semi-Discrete Central-Upwind Schemes for Hyper- 
bolic Conservation Laws and Hamilton- Jacobi Equations , SIAM J. Sci. Comp., 23 , (2001), 
pp. 707-740. 

[10] Levy, D., Puppo, G., Russo, G., A Fourth Order Central WENO Scheme for Multi- 
Dimensional Hyperbolic Systems of Conservation Laws , SIAM J. Sci. Comp., 24 , (2002), 
pp. 480-506. 



Modeling Dynamics of the Solar Atmosphere 


31 


[11] Lites, B. W., Skumanich, A., A Model of a Sunspot Chromosphere Based on OSO 8 Obser- 
vations , ApJS, 49 , (1989), pp.293-316. 

[12] Liu, X.-D. L., Osher, S., Nonoscillatory high order accurate self-similar maximum principle 
satisfying shock capturing schemes I , SIAM J. Numer. Anal., 33, no. 2 (1996), pp. 760-779. 

[13] Shibata, K., Suematsu, Y., Why are Spicules Absend Over Plages and Long Under Coronal 
Holes?, Solar Physics, 78 , (1982), pp. 333-345. 

[14] Sterling, A. C., Solar Spicules: A Review of Recent Models and Targets for Future Obser- 
vations, Solar Phys., 196 , (2000), pp. 79-111. 

[15] Sterling, A. C., Mariska, J. T., Numerical Simulations of the Rebound Shock Model for Solar 
Spicules, ApJ, 349 , (1990), pp. 647-655. 

[16] Suematsu, Y., Shibata, K., Nishikawa, T., Kitai, R., Numerical Hydrodynamics of the Jet 
Phenomena in the Solar Atmosphere I. Spicules, Solar Physics, 75, (1982), pp. 99-118. 

[17] Sweby P.K., High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation 
Laws, SIAM J. Numer. Anal., 21 , (1984), pp. 995-1011. 

[18] Vernazza, J. E., Avrett, E. H., Loeser, R., Structure of the solar chromosphere. Ill - Models 
of the EUV brightness components of the quiet-sun, ApJS 45, (1981), pp. 635-725. 


